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ABSTRACT 

We study the effect of the non-linear process of ambipolar diffusion (joint transport 
of magnetic flux and charged particles relative to neutral particles) on the long-term 
behavior of a non-uniform magnetic field in a one-dimensional geometry. Our main 
focus is the dissipation of magnetic energy inside neutron stars (particularly magne- 
tars), but our results have a wider application, particularly to the interstellar medium 
and the loss of magnetic flux from collapsing molecular cloud cores. Our system is a 
weakly ionized plasma in which neutral and charged particles can be converted into 
each other through nuclear beta decays (or ionization-recombination processes). In 
the "weak-coupling" limit of infrequent inter-particle interactions, the evolution of 
the magnetic field is controlled by the beta decay rate and can be described by a non- 
linear partial integro-differential equation. In the opposite, "strong-coupling" regime, 
the evolution is controlled by the inter-particle collisions and can be modelled through 
a non-linear diffusion equation. We show numerically that, in both regimes, ambipolar 
diffusion tends to spread out the magnetic flux, but, contrary to the normal Ohmic 
diffusion, it produces sharp magnetic field gradients with associated current sheets 
around those regions where the magnetic field is weak. 

Key words: diffusion - plasmas - magnetohydrodynamics (MHD) - ISM: magnetic 
fields - stars: magnetic fields - stars: neutron 



1 INTRODUCTION 

Ambipolar diffusion is the joint drift of charged particles and 
the associated magnetic flux with r espect to the neutral par - 
ticles in a partially ionized plasma. iMestel fc Spitzen (|l956l ) 
first proposed it in order to explain the loss of magnetic flux 
from the dense cores of molecular clouds, required for the 
formation of stars, starting an active field of research in this 
area. Later, it was suggested to also pl ay a role in the decay 
of the magnetic fields of neutron stars |joned|l98^ : iHarrisonl 
lf99ll : |Petiiickll992l : lGoldreich fc ReiseneggeJl99a , hereafter 
GR-92) which became particularly relevant with the iden- 
tification of "magnetars", neutron stars whose main power 
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In a previous paper 

l2008l : hereafter Paper I), 



we established a multifluid for- 
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malism in which it is possible to study the lo ng-term evo- 
lution of magnetic fields in neutron stars (see iReiseneggeij 
I2OO9I for a discussion of the main properties of the magnetic 
field equilibria and their subsequent long-term evolution). 
In that work, and following the ideas developed by GR-92, 
we included the effects of several physical processes that are 
also relevant for star formation and protoplanetary disks, 
including ambipolar diffusion, Hall drift (non-dissipative ad- 
vection of the magnetic field by the associated electrical cur- 
rent), and ohmic diffusion (dissipation of currents through 
the electrical resistivity). 

Here we continue this study and concentrate our anal- 
ysis on the long-term evolution of the magnetic field 
caused by ambipolar diffusion aided by beta decays. Fol- 
lowing the same philosophy of Paper I, and as a first 
approach to the understanding of our general formal- 
ism, we focus on a simplified, one-dimensional configu- 
ration in which the magnetic field points in one Carte- 
sian direction z but varies only along an orthogonal di- 
rection X, i. e., B — B(x,t)z. Such models have also 
been considered in several studies of ambipolar diffu- 
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sion i n the interstellar medium (iMouschovias fc Paleologoul 
Il98ll : [siiul 1 19831 : [Brandenburg fc Zweibej 1 19941 ). although 
some of the assumptions differ from case to case. In our 
analysis, we consider separately two relevant limits, sim- 
ilar in spirit, though not exact ly equivalent to those of 
IMouschovias fc Paleoloeod (|l98ll l: 

• In the weak coupling limit, there are few collisions be- 
tween the particles, and the beta decays proceed slowly. 
Therefore, the particles can reach the diffusive equilibrium 
easily, but it takes much longer to reach the chemical equi- 
librium. 

• In the strong coupling limit, there are many collisions 
between the particles, and the beta decays proceed fast, 
so the (local) chemical equilibrium is reached much more 
quickly than the diffusive equilibrium. 

For each of these cases, we find that the long-term evo- 
lution of the magnetic field can be modelled by a single 
equation that gives the time-derivative dB/dt at a given 
instant t only in terms of the configuration of the magnetic 
field at the same instant, B{x, t). This makes it easy to carry 
out numerical simulations of the evolution of some selected 
non-linear magnetic field profiles, and even find some exact, 
analytical solutions. 

In Sect. [2] of this paper, we briefiy review the one- 
dimensional model of neutron star magnetic field evolution 
introduced in Paper I, paying particular attention to its 
characteristic evolutionary time-scales. 

In Sects. [3] and [l] we obtain the equations for the long- 
term, asymptotic magnetic field evolution promoted by am- 
bipolar diffusion in each of the two opposite regimes men- 
tioned above, and we make numerical simulations of the evo- 
lution of different initial magnetic field configurations. We 
show that, in both cases, the magnetic flux of a given sign 
tends to spread out, but singularities develop at the null 
points w here regions with different signs meet, as previously 
found bv iBrandenburg fc Zweibel (|l994 ). In the weak cou- 
pling case, these singularities correspond to current sheets 
that are dissipated by resistive effects, in this way leading 
to reconnection. In the strong coupling case, the singular- 
ities have a somewhat different character (a smoothly di- 
verging current density) and might lead directly to recon- 
nection even in the case of no ohmic resistivity (but see 
iHeitsch fc Zweibelll2003al lbh. Finally, in Sect.E we give the 
main conclusions of our study. 



2 ONE-DIMENSIONAL MODEL AND 
EVOLUTIONARY TIME-SCALES IN 
NEUTRON STARS 

We model the neutron star interior as an electrically neu- 
tral and slightly ionized plasma composed of three mobile, 
strongly degenerate, particle species: neutrons (n), protons 
(p), and electrons (e). We account for binary collisions and 
weak interactions (causing nuclear beta decays) between the 
particles, and allow for strong interactions between neu- 
trons and protons by writing each of their chemical po- 
tentials as a function of both of their number densities: 
fJ"n,p = ^'n,p{nn,np), wMle considering the electrons as an 
ideal, relativistic Fermi gas, whose chemical potential is a 
function only of their number density, ~ Heine). 



We study a one-dimensional geometry in which the 
magnetic field points in one Cartesian direction z, but varies 
only along an orthogonal direction x as B{f,t) = B{x,t)z, 
and assume that all physical variables vary only along x. 
Since, in neutron star conditions, the ratio of the magnetic 
pressure B'^/Stv to the pressure of the charged particles is 
very small, the magnetic force causes only small pertur- 
bations to the hydrostatic equilibrium state of the non- 
magnetized star. 

These assumptions are generally not true in molecular 
cloud cores, where the ionization fraction tends to be ex- 
tremely low, whereas the magnetic field can be near equipar- 
tition with the neutral gas pressure. In this sense, our deriva- 
tion will be valid only for the case of neutron stars, although 
we will see that some of the results agree with those of other 
authors, obtained under somewhat different assumptions. A 
more general treatment appears to be difficult and not to 
yield simple results. 

For the reasons stated, we consider a non-magnetized, 
fixed background system in hydrostatic and chemical equi- 
librium and introduce small perturbations to the number 
density of each species i as ni{x,t) = noi{x) + Sni{x,t), with 
the subscript zero labelling the background number densities 
and \5ni{x,t)\ <^ nQi{x). The associated chemical potential 
perturbations are given by /li = fioi + Sjj,i. 

The long-term magnetic field evolution implies small 
particle velocities that change over long time-scales, much 
longer than the very short dynamical times that are only 
relevant shortly after the formation of the star (i.e, sound or 
Alfven time-scales). Therefore, at all times we use a slow- 
motion approximation in which we neglect the acceleration 
terms in the equations of motion for the particlesQ 

Taking account of all these considerations, the system 
of non-linear partial differential equations governing the evo- 
lution is (see Paper I for the derivation) 
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^ This approximation should be roughly true also for 
star f ormation, and was indeed assumed by some au- 
thors jsiiul Il983l: iBrandenbure fc Zweibej Il994l l, but not by 
iMouschoviasfcPaieologoiJ l|l98lh . who dropped the charged- 
particle pressure, but kept the inertial terms, and therefore ob- 
tained a differential equation of second rather than first order in 
time for the magnetic field evolution. 
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The magnetic field evolution is given by Eq. ([T]), where 
Vc denotes the velocity of the charged particles (the same 
for electrons and protons), Vn is the velocity of the neu- 
trons, and va = Vc — Vn is the ambipolar diffusion veloc- 
ity0 The first term on the right-hand side of this equa- 
tion causes an advection of the magnetic flux with a ve- 
locity Vc, while the second term describes the ohmic dif- 
fusion of the field, where ctq is the electrical conductivity. 
In neutron star core conditions, the electrical conductivity 
is (70 ~ 10^*s^^, thus, the evolution of the large-scale mag- 
netic field through ohmic diffusion proceeds very slowly, with 
a time-sc ale tphmic ^ 10^^ yr, longer than the age of the 
universe (|Bavm. Pethick. fc Pin"e3ll969l '). Qualitatively sim- 
ilar conditions hold in essentially all astrophysical settings. 
Therefore, in the rest of this paper, we neglect the ohmic 
term in the magnetic field evolution equation and focus only 
on the magnetic evolution due to the advective term. How- 
ever, the time-scale for ohmic diffusion scales with the square 
of the characteristic length of the magnetic field variations 
(see Paper I), thus, the influence of ohmic dissipation can 
be important in regions with strong spatial magnetic field 
variations. 

The evolution of the particle number density pertur- 
bations, for the different species, is given by Eqs. ^ and 
((31), where Sric — Sup = Suc is the perturbation of the 
charged particle number density keeping charge neutrality 
and Sub = i5w„ -I- Suc is the perturbation of the baryon 
number density. Weak interactions cause beta decays (con- 
version of charged particles into neutrons and the opposite) 
that tend to reduce deviations from the chemical equilib- 
rium state between charged particles and neutrons, with a 
net rate coefficient A and Sfic = Sfic + Sfip. The chemical 
equilibrium is achieved when Sjic = Sfi^. 

Continuing our description, Eq. ((5]) shows that the 
ambipolar diffusion velocity va is controlled by the col- 
lision rate between the charged particles and neutrons, 
which is proportional to the parametre yen- This veloc- 
ity is driven by the Lorentz force but c hoked by the 
char ged particle pressure gradient it produces (jPethickl 19921 : 
iGold rcich & RciscncsEcr 1993). 

Regarding the neutrons, we see from the right-hand side 
of Eq. Q that they move with a velocity Vn as long as 
the system is not in a state of magneto-hydrostatic equi- 
librium. In this state, the Lorentz force is balanced by the 
pressure gradients of all the particles (charged particles and 
neutrons). We also see that the neutron velocity is controlled 
by the parametre q, whose meaning we explain in the next 
paragraph. 

Since we are interested in a numerical solution of our 
equations, we have to take into account that modelling the 
true dynamical time-scales (sound or Alfven time-scales of 
milliseconds to seconds) would require a time step many or- 
ders of magnitude shorter than that required to simulate the 
long-term evolution, making the simulation computationally 
unfeasible. For this reason, in Paper I we introduced a slow- 



^ Strictly speaking, all these are the x-components of the respec- 
tive velocities, which cause changes in the distributions of parti- 
cles and magnetic flux. The charged-particle velocities must also 
have a i/-component, responsible for the current that acts as the 
source of the magnetic field. 



motion approximation, neglecting the acceleration terms and 
instead introducing a small, artificial, friction-like force act- 
ing on the neutrons (the most abundant species) of the form 
—no„av„, where the parametre a is chosen in such a way 
that the time to reach magneto-hydrostatic equilibrium is 
long enough for the numerical code to be able to deal with it 
(and therefore, much longer than the real dynamical times), 
but shorter than the long time-scales over which the mag- 
netic field evolves. We showed in Paper I that the latter 
time-scales, which are the astrophysically interesting ones, 
are unaffected by the choice of a. 

We require the conservation of both the magnetic fiux 
$ = Jq B dx and the baryon number perturbation 5Nb = 
Jq Sub dx during the evolution of our system, which spans 
the segment ^ x ^ d. We ensure this through the bound- 
ary conditions 

Vcix = 0,t)=Vcix = d,t) = 0, (6) 

Vnix ^ 0,t) = Vnix = d,t) ^ 0, (7) 



^ix^O,t)^-ix = d,t) = 0. (8) 

The system of Eqs. (|H3|) describes the evolution of three 
coupled variables: the magnetic field, the charged particle 
density perturbation, and the baryon density perturbation. 
In Paper I, we estimated the three associated characteristic 
evolutionary time-scales of this set of equations correspond- 
ing to exponentially decaying eigenmodes in the linear ap- 
proximation and showed that they characterize the approach 
to three successive equilibrium states. Here we briefiy de- 
scribe the evolutionary stages and summarize the relevant 
time-scales with the main goal of establishing the basic ideas 
that will be used in our subsequent analysis (see Paper I for 
more details about their derivation). 

The shortest time corresponds to the approach of the 
magneto-hydrostatic equilibrium, controlled by a in our 
model, as already explained. 

In a longer time, two alternative processes compete: 

(i) The particle species move relative to each other, con- 
trolled by the inter-particle collisions, in a tendency to reach 
a diffusive equilibrium state in which the fluid forces acting 
on each species separately are balanced. 

(ii) The weak interactions convert particles from one 
species into another, tending towards a chemical equilibrium 
state. 

If the relative motion of charged and neutral particles 
proceeds much faster than the conversion from one into an- 
other, a diffusive equilibrium is achieved in the system, char- 
acterized by the balance equations 

and 

= nOn^On^ ( I . (10) 

dx \^J.On ) 

The inter-particle collision frequency controls the time-scale 
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to reach the diffusive equihbrium, which is given by Paper 
I, 



idrag 



V SnoB J , 



4.5 X 10"' Ll r| yr. (11) 



where Lg = i/(10^cm) and Ts = T/(10**K). If the beta 
decays are much faster, the chemical equihbrium state 5/i„ = 
Sfic is reached in a time-scale controlled by the beta decay 
rate, 



4.3 X 10^ Tg'^ yr. (12) 



Thus, we have two relevant limits: 

(i) In the weak coupling limit, tdrag ^ tweak and the weak 
interactions operate slowly, therefore, the system relaxes to 
a diffusive equilibrium in a time-scale tdrag, but remains out 
of chemical equilibrium. 

(ii) In the opposite, strong coupling limit, tdrag ^ tweak, 
the inter-particle collisions are very frequent which impedes 
a fast achievement of the diffusive equilibrium, while the 
system can relax to the chemical equilibrium in the time- 
scale i^jjeafc- 

The transition between these limits is achieved when tdrag ~ 
tweak- From Eqs. pip and p2|) we infer that the condition 
for this transition is Ts ~ 5.6 Lj which gives tdrag ~ 



til 



14 L 



3/2 



yr- 



Since the density perturbations are assumed to be small, 
these will not involve large motions of the particles, and 
therefore not cause a substantial change in the magnetic 
flux distribution. Note, however, that these two equilibria 
are incompatible with each other as long as a spatially non- 
uniform magnetic field is present, therefore a full equilibrium 
will only be reached in a much longer time-scale, on which 
the magnetic field is made uniform (in our model) or ex- 
pelled from the system (likely more realistic in a true astro- 
physical setting). For the determination of this much longer 
time-scale, on which the magnetic field does evolve substan- 
tially, we again consider the two opposite regimes discussed 
above. In the weak coupling limit, the system reached the dif- 
fusive equilibrium but not the chemical equilibrium, during 
the previous stage. During this much longer stage, the weak 
interactions slowly convert charged particles into neutrons 
in a tendency to reduce the charged-particle pressure gra- 
dient that counterbalances the magnetic pressure gradient. 
This causes a slight deviation from the diffusive equilibrium, 
producing a joint transport of the charged particles and the 
magnetic flux at a small ambipolar diffusion velocity va, 
always keeping the system very close to diffusive equilib- 
rium. This interplay continues until both the pressure and 
magnetic fleld gradients dissapear, which occurs in a long 
time-scale that depends on the magnetic field strength and 
the weak interaction rate: 



(weak) 
ambip 



Pt^ 



AB2 



1.7 X 10^ Tg-^ yr, (13) 



with /3 = SvrnocJiOn (9/ion/9nos)„^^ /-B^ » 1, which is 
roughly the ratio of the charged particle pressure to the 
magnetic pressure and we defined Bis = 5/(10^^^0). 

In the strong coupling limit, the relative motion of the 



charged particles and neutrons is strongly suppressed by the 
inter-particle collisions, which delay the diffusive equilib- 
rium state, while in comparison the chemical equilibrium 
is reached quickly. The deviation from diffusive equilibrium 
promotes, as before, a joint motion of the charged particles 
and the magnetic flux, with a very small ambipolar diffusion 
velocity va- This movement yields the diffusive equilibrium 
in a very Ion time-scale controlled by the collision rate be- 
tween particles and by the magnetic-field amplitude. 



{drag) 
ambip 



tdrag 



~ 1.8 X 10=^ B^^ Ll Ti yr. 



(14) 



Note that at the transition between the strong and 
weak coupling limits (Tg ~ 5.6 L^^^''), the ambipolar dif- 
fusion timescales are of the same order: t^^^f^ ~ t^ZnUp ~ 
5.5 X 10'' B^l^ yr. This value corresponds to the short- 
est possible ambipolar diffusion time, since the relevant 
timescales increase both towards higher and lower temper- 
atures. Although the Hall effect is not present in our one- 
dimensional calculations, it is important to asses its poten- 
tial importance relative to ambipolar diffusion. Consider its 
timescale,: 



tuall 



4nnoeL^ 
cB 



3 X 10" B. 



-1 r 2 

5 yr. 



(15) 



compared to the minimum ambipolar diffusion 

— 1/2 

timescale estimated above. We find that, if Bis < 0.2 , 
the Hall drift is likely to play a dominant role. It might not 
be important in magnetars, where Bis > 1, except possibly 
in reconnection layers, where is «C 1- 

In the present paper, we are interested in the details of 
the evolution of the magnetic field in the strong and weak 
coupling limits. Therefore, in the next sections we will ob- 
tain the differential equations modelling the evolution of the 
magnetic field in the time-scales given by Eqs. (|13[) and {2)). 



3 MAGNETIC FIELD EVOLUTION IN THE 
WEAK COUPLING LIMIT 

3.1 Derivation of the asymptotic evolutionary 
equation 

In this section, we derive a single integro-differential equa- 
tion that models the magnetic field evolution in the weak 
coupling limit {tdrag ^ t^^ak), in which neutral and charged 
particles drift easily with respect to each other, and the main 
bottleneck is the (slow) rate at which they can be converted 
into each other, in order to decrease the charged-particle 
pressure gradients that balance the Lorentz force, imped- 
ing the magnetic fiux to spread. This limit is not likely to 
be relevant in the interstellar medium (|Shull 19831 ) ■ although 
it roughly corresponds to one of the limits considered by 
iMouschovias fc Paleologoul l|l98ll ). In neutron stars, it be- 
comes important at later stages, once their temperature is 
low enough. 

Hereafter, we take the background properties to be ho- 
mogeneous, i.e. the variables with a subscript do not de- 
pend on the position. We assume that the system has al- 
ready reached the diffusive equilibrium, which also implies 
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the magneto-hydrostatic equiUbrium. However, since weak 
interactions are slow, the system may be out of chemical 
equilibrium, namely, S^n 7^ S/^c- From Eq. (|10p we infer a 
spatially uniform distribution of the neutrons, 



S/in = g{t), 



(16) 



with g{t) an arbitrary function that depends only on the 
time variable. From Eq. ((9|, the diffusive equilibrium of the 
charged particles implies 



Stt 



nocSfic + — = fit)- (17) 
Since ISricl <C noc and the time-scale for the evolution 



is t 



(weak) 



L/vc, we can compare the terms in Eq. 



dUr 



dt 



SUc 



,(weak) 
ambip 



nocVc 



d 



(nocVc) 



(18) 



so we can neglect the time derivative and write this equation 
as 



dvc A 

— — = (dfic — O/ln 

ox noc 



(19) 



Integrating this equation between a; = and an arbitrary 
internal point with coordinate x, we obtain 



Vc{x,t) 



X 



[Sfj.c{x' ,t) — Sfin{x' ,t)] dx' , (20) 



where we used the boundary condition Vc{0, t) — 0. Replac- 
ing Eqs. ([T6|) and 1(17} in Eq. (f20|) and using the boundary 
condition Vc{d,t) = to eliminate f{t) and g{t), we obtain 



Vc{x,t) 



X 



x{d—x) 



dx' 



8^" 



dx' 
(21) 

Note that the velocity at any given point is proportional to 
the parametre X controlling the weak interaction rate and 
to the difference of the average magnetic pressure to the left 
and to the right of this point. This is because, in diffusive 
equilibrium, a stronger magnetic pressure corresponds to a 
lower density of charged particles, and the magnetic flux 
can only spread as the reactions modify the charged-particle 
density. 

Replacing Vc in Eq. ([T]), we obtain the integro- 
differential equation governing the magnetic field evolution 



dn^^ dx 









x{c 


i-x) 


( ' f 

\d-x 







El 

Stt' 



dx' 



8^' 



dx' B 



(22) 



The characteristic time-scale from this equation is the 
same as estimated in Paper I. If we scale the position variable 
x in Eq. (|22p to the total length of the system d, the time 



variable to t^^up, and the magnetic field to some character- 
istic magnetic field (for our subsequent numerical analysis it 
is chosen as the maximum of the initial magnetic field pro- 
file) , we can write this equation with dimensionless variables 
as 



dB 
'dt 



d_ 

dx 



x{l 



B^dx' 



B^dx' B 



(23) 



3.2 Analytical and numerical solutions 



Below, we solve Eq. ((23} by using a immerical finite- 
difference scheme. However, and with the goal of testing 
the adequate performance of the numerical method, it is de- 
sirable to know some analytical solutions to compare with 
the numerical results. Following this philosophy, we try an 
analytical solution of Eq. H23|) in the form of a "step" pro- 
file, motivated by the fact that we observed it as a generic 
asymptotic state in the evolution of some of the initial pro- 
files that we study in the following paragraphs. Thus, we 
make the Ansatz 



B{x,t) = 





Bs{t) 



<x < xi{t), 

Xl{t) < X <1. 



(24) 



Using the condition of flux conservation, <I> = Bs{t){l — 
xi{t)) = constant, and replacing Eq. (|24|l in the dimension- 
less version of Eq. ((21} evaluated at x = x\, we get 



dx\ 
~dt 



X\ 



Integrating Eq. (|25|) yields 



t - to = 



1 

" $2 



In 



Xl{t) 

Xl{tQ) 



{xi{t) — Xl{to)) 



(25) 



(26) 



where to is some reference time. Inverting this equation gives 
xi{t), while Bs{t) can be obtained through the condition of 
flux conservation. 

The step profile Eq. (|24p presents a sharp gradient 
around xi{t). Since we are neglecting the ohmic dissipation, 
any attempt that we made to use this step profile as an 
initial condition in the finite-difference scheme results in a 
numerical instability. In order to overcome this difficulty, we 
compare in Fig. (T] the analytical evolution of Eq. (|24l) with 
the finite-difference solution of Eq. (|23p . using as initial con- 
dition a smooth profile. 



B{x,0) 



1 + tanh 



X — Xo 



(27) 



which is a good approximation to Eq. (|24|) . 

For this comparison, the initial parametres xi{to = 0) 
and Bs{t = 0) are chosen so that the step profile given by 
Eq. (|24|l and the smooth profile given by Eq. (|27|) share 
the same initial magnetic flux. We see from Fig. [T] that, 
except for the differences around xi{t), the evolution of both 
profiles is similar. 

Now, we intend to verify whether the asymptotic evolu- 
tionary equation Eq. (|23|) correctly describes the long-term 
behavior of the magnetic field. This requires a comparison 
of the asymptotic evolution of the magnetic field as given by 
the full set of Eqs. ((T}, ((2]), and ((3]) with the corresponding 
evolution that arises when using only Eq. ((23}. In Fig. (2] we 
carry out this comparison for the evolution of a Gaussian 
initial magnetic field profile given by 



B{x,0) = exp[- 



-S (X - 



- Xo 



(28) 



As initial condition for the particle densities we set 
SnB{x,0) = (5nc(a;,0) = 0. Note that this comparison is 
valid at late times, when this profile has evolved to be con- 
sistent with the diffusive equilibrium. So, we start this com- 
parison from the instant labelled with number (2) in this 
figure, which corresponds to t* = 2^tdrag- We see that the 
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Figure 1. The analytic evolution of an initial magnetic step pro- 
file of the form given by Eq. I|24| l (thin line), and comparison with 
the finite-difference evolution of an initial profile given by Eq. i27l 
but using Eq. I|23| l (thick line). The time scaling in this figure and 



Aweak) 



ones in this regime is such that , 



1 and 



in the following 

the time progression is labelled as: (1) t = 0, (2) t = 2.45, (3) 
t = 8.07. We used the parametres xq = 0.6, to = 0, a = 0.04. 
Bs{0) = 1 and xi(0) = 0.6 were chosen so that both profiles 
share the same magnetic flux. 
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Figure 2. The evolution of a Gaussian initial magnetic field pro- 
file given by Eq. I I28I I with s = 20, and xg = 0.5 from the finite- 
difference solution of the coupled system of Eqs. JTJ, l[2]l, and Q 
(thick line). In the instant t = 20tdrag [labelled with (2)] this pro- 
file has evolved to be consistent with the diffusive equilibrium. So, 
from this instant we also follow the evolution given by the numeri- 
cal finite-difference solution of Eq. I I23I I (line with empty squares). 
In addition, from the instant t = 50t'™\'''^^ [labelled with (3)1 , we 

' amoip L \^ /J ' 

also compare with the box-type solution given by Eq. I I29I1 (thin 
line), where xi at this instant was calculated so that all the pro- 
files share the same magnetic flux, using the flux conservation 
relation $ = -83(1 — 2xi). For this simulation, we set as initial 
condition for the particle densities 6nB{x,0) = Snc{x,0) = 
and the parametres: noc/no„ = 0.04, L/d = 0.08, /9 = 2.0, and 
t-weak/tdrag ~ 100. The labels in the figure represent the instants: 
(1) t = 0, (2) t = 20taraa = 0.1, (3) t 
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100. 



Figure 3. The evolution of Eq. I I32II from the finite-difference 
solution of the full coupled system of Eqs. Q, ((2]|, and ^ (full 
line). Similarly to Fig. (2] in the instant t = lOt^rag [labelled 
with (2)], this profile has evolved to be consistent with the diffu- 
sive equilibrium. So, from this instant we also follow the evolu- 
tion given by the numerical finite-difference solution of Eq. l|23p 
(points). For this simulation we set again as initial condition 
5ng{x,0) = 5nc{x,0) = and the same background parame- 
tres of Fig. [2] We set (3 = 1.1 and t^^ak/tdrag ~ 100. The labels 
in the figure represent (1) t = 0, (2) t = lOt^rag = 0.09072, (3) 
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asymptotic behavior of the magnetic field (labels (3) and 
(4) in the figure) calculated from these two methods is the 
same. With the purpose of characterizing this asymptotic 
behavior analytically, we construct an analytic "box" -type 
solution of the form 



B{x,t) = 




< a; < xi{t), 

Xi{t) < X < 1 — Xi{t), 

1 - Xl{t) < X <1. 



(29) 



with 



dxi 
~dt 



Xl 



2xi 



where the magnetic flux is "1> = Bs(t)(l 
is given by 



(30) 

2xi{t)), and t{xi) 



$2 



In 



Xljt) 
Xl(to) 



- 2{xi{t) - xi{to)) 



(31) 



For the instants (3) and (4) in Fig. [2] we show a full line 
that represents the solution given by Eq. ((29}. We see that 
this analytical characterization works well. 

In Fig. [3] we show the evolution of a harmonic initial 
magnetic field profile 



B{x,0) = — cos(7ra;) 



(32) 



In this figure, we compare the evolution of this pro- 
file given by Eq. (|23p with that given by the fuU system 
of coupled Eqs. ((2]), and For the latter, the ini- 
tial conditions for the particle density perturbations are 
5nB{x, 0) — Snc{x, 0) = 0. At the instant labelled with num- 
ber (2), t* — lot drag, tMs initial profile has evolved to be 
consistent with the diffusive equilibrium, therefore, from this 
instant we also calculate the evolution of the magnetic field 
through Eq. (|23p . We see in this figure that the asymptotic 
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Figure 4. Evolution of a double Gaussian profile of the form 
B(x,0) = exp[-400(x - 0.2)2] _,_ q 7 

exp[— 100(a; — 0.8)'^], using 
Eq. II23I I. The labels indicate the time progression for different 
instants, (f ) t = 0, (2) i = 4.0, (3) t = 60.0, (4) t = 400.0. 
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Figure 5. Evolution given by Ea.l l23l l for the same instants as in 
Fig.U but for an initial profile with one of the Gaussians inverted: 
B(a:,0) = exp[-400(a; - 0.2)2] _ 0.7 exp[-100{a; - 0.8)^] 



behavior of the magnetic field [instant (4)] calculated from 
these two methods is the same, which again verifies the va- 
lidity of Eq. H23|) to model the asymptotic evolution of the 
magnetic field. 

In summary, the last results confirm the adequacy of 
Eq. (|23p to describe the asymptotic evolution of the mag- 
netic field promoted by ambipolar diffusion but controlled 
by beta decays. So, in what follows we explore the evolu- 
tion of different magnetic profiles using only this equation. 
This also has the advantage that we can increase the time- 
step with respect to that needed when solving the full set 
of Eqs. (HI, (d}, and ((3)l without generating numerical insta- 
bilities. In Figs. we show the evolution of different ini- 
tial profiles. From the numerical results, we can infer some 
generic properties: 

• The ambipolar diffusion process tends to spread out the 
magnetic fiux and the total magnetic fiux is conserved as it 
is expected from our boundary conditions. 

• Contrary to the normal Ohmic diffusion, in Fig. [S] it 
is observed that the magnetic field is not smoothed over 
all the space. It becomes uniform only over regions whose 
boundaries are points where the magnetic field is zero, and 
across which it jumps between values of the same magnitude, 
but opposite signs, making the magnetic pressure jif^-ii) 
equal on both sides. 

• The magnetic field nulls move according to the local 
values of the charged particle velocity, as given by Eq. (|2ip . 
This velocity is a continuous function of x, so the flux can 
be spread out or compressed, but no reconnection (mutual 
elimination of opposite field lines) occurs in the absence of 
Ohmic dissipation, contrary to the strong coupling limit to 
be studied in the next section. Of course, the formation of 
steep gradients makes it possible for Ohmic dissipation and 
therefore reconnection to occur in a realistic astrophysical 
setting. 



4 MAGNETIC FIELD EVOLUTION IN THE 
STRONG COUPLING LIMIT 

4.1 Derivation of the asymptotic evolutionary 
equation 

In this section, we study the magnetic field evolution in the 
opposite limit of the last section, namely, the strong coupling 
limit, tdrag 3> t^jeak, in which the conversion from charged to 
neutral particles and vice- versa is essentially instantaneous, 
but their relative motion is impeded by a strong mutual 
coUisional drag force. This limit is relevant for young, hot 
neutron stars, as well as in the in terstellar medium (|Shul 
1 19831 : [Brandenburg fc Zweibellll994l ) . 

In order to study the long-term evolution, we con- 
sider that the system has reached both the magneto- 
hydrostatic and chemical equilibria. From Eq. Q, the 
magneto-hydrostatic equilibrium condition implies 



B 

no„Sfi„ + nocSfic + 7;— ^ h{t), 
on 



(33) 



where h{t) is an arbitrary function depending only on the 
time variable. On the other hand, the condition for chemical 
equilibrium implies 



5fJ,n — 5fJ,c = 

Combining Eqs. (|33p and (|34|l . we obtain 

h(t) 



= S^n 



(34) 



(35) 



As done in Eq. (|18p . we can compare the time derivative to 
one of the spatial derivatives in Eq. ((2]) , 



dSn£ 



dt 



Sub ^ npcVc 



.(drag) 
amhip 



noc 



dvc 



dx 



(36) 



Thus, we can neglect the time derivative of Sub in compar- 
ison with the spatial derivative term in Eq. ((2]) and write 



non— — = —noc 



Integrating Eq. (I37|) between x — and an inner point 
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X with the boundary conditions Vc{Q,t) = ii„(0,t) = 0, and 
using Eq. (|35p . we obtain 



d 



(38) 



We see from Eq. p8|) that the Lorentz force drives the mo- 
tion of magnetic flux and charged particles at the velocity 
Vc, which is controlled by the inter-particle collisions through 
the factor After including Eq. (|38p in Eq. we ob- 

tain the equation for the magnetic field evolution as 



dB 
'dt 



d_ 

dx 



(vcB) 



non 



d 



8nnocnlg-fcn dx 



B 



dB^ 
dx 



(39) 



Scaling the time variable to t^'^^j^' and the position vari- 
able to the total length of the system d, this equation can 
be written in dimensionless form as 



dB_ 

at 



d_ 

dx 



B 



dB^ 

dx 



(40) 



Eq. (|40p belongs to a group of non-linear diffusion equa- 
tions or porous medium equations whose mathematical prop- 
erties have been stud i ed by differ ent authors (see e.g. iTuckl 
1 19761 : IVasauedlioOTl ). Isiiul (|l983l ) obtained it for ambipo- 
lar diffusion in interstellar gas, in a similar, strong-coupling 
regime as considered here, but ignoring the pressure of the 
charged particles, while allowing for arbitrarily large den- 
sity perturbations, and using a Lagrangian coordinate sys- 
tem moving with the neutral fluid. It was re-derived under 
slightl y different assumptions by [Brandenburg fc Zweibell 
l| 19941 ). 



4.2 Analytical and numerical solutions 

4-2.1 Exact analytical solution 



As an analytical solution of Eq. (|40[) tha t will allow us to 
validate our numerical results, we foUow iTuckl (|l976l ) and 
consider 



B{x,t) 







Bs(t)Jl- 



a,(t) 



if 
if 



\x — Xo\ > as{t) 
\x — xo\ ^ as{t). 



(41) 

This solution is a semi-ellipse in the {x, B) plane, cen- 
tered at (xojO), with semi-axes as(t) and Bs{t). Replacing 
Eq. (|41[) in Eq. (|40p and using the flux conservation condi- 
tion $ = (ji /2)Bs{t)as{t) = constant, we obtain the differ- 
ential equation for as{t), namely. 



dus 



9 ^ ' 



with the solution 



■ 8B?(£) 



-,1/4 



(42) 



(43) 



where t* is a reference time. 

In Fig. |6l we compare this analytical solution to the 
corresponding evolution obtained from the finite-difference 
numerical integration of Eq. (|40p . 

We see that there is a good agreement between the two 
methods except at the last instant, labelled with (4) in the 
figure. This is expected, since at this instant the elliptic 
part of the profile has crossed the boundaries of the system, 
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Figure 6. Evolution of an initial "elliptic" profile of the form 
Eq. l|41|l from Eqs. Ij42p . l|43|l (line), and evolution given by tlie 
numerical finite-difference solution of Eq. II40I I (empty squares). 
The time instants are (1) t = 0, (2) i = 0.02, (3) t = 0.3, (4) 
t = 0.9 



so the analytical solution is no longer valid. Note also that 
the latter does not conserve magnetic flux once the ends 
of the ellipse have reached the boundaries. We also see at 
the instant labelled with (4) how the magnetic fleld given by 
Eq. H40[) is almost homogeneously distributed across the sys- 
tem, as expected from the magnetic-flux conserving bound- 
ary conditions that we are using for this equation. Thus, the 
flnite-difference method that we are using to solve Eq. (1401) 
appears to be working well. 



4-2.2 Numerical explorations 

Hereafter, and with the main goal of establishing a compar- 
ison, we follow the evolution of the same initial proflles that 
we used in the previous section, where we dealt with the 
opposite regime. 

In Fig. [T] we show the evolution of an initial Gaussian 
magnetic field profile given by Eq. (|28|) . as given by the full 
coupled system of Eqs. ((T]), and ^ with 5ns (x,0) = 
5nc{x,0) = 0. At the instant (2), t* — lOt^eak, the system 
should be very close to chemical equilibrium. Thus, from 
this instant onwards, Eq. (|40p should be valid, so we also 
follow the evolution of the magnetic field obtained from the 
finite-difference solution of this equation. We see that the 
asymptotic behavior of the magnetic field (labels (2) and 
(3) in the figure) calculated from these two models is the 
same. In addition, for the instants (2) and (3), we show for 
comparison the analytical solution [Eq. (|41|) ]. 

We see that this explicit solution is an adequate descrip- 
tion of the asymptotic evolution. 

The previous results confirm the adequacy of Eq. ((40)) 
in describing the asymptotic evolution of the magnetic field 
promoted by ambipolar diffusion and controlled by inter- 
particle collisions. We now examine the evolution of different 
magnetic profiles using only this equation. In Figs. 1 81101 we 
see the evolution of different magnetic field profiles. From 
these numerical results, we can infer the following properties: 

• As in the weak coupling limit, the ambipolar diffusion 
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Figure 7. Evolution of a Gaussian initial field profile of the form 
given by Eq. I|28|l with s = 6, and xq = 0.5, from the finite- 
difference solution of the coupled system of Eqs. Q, and Q 
(line with points). For this simulation we set as intitial condition 
for the particle densities 5ng{x, 0) = <5nc(x, 0) = 0. In the instant 
t = lOtiueafe [labeled with (2)], this profile has evolved to be con- 
sistent with the chemical equilibrium. So, from this instant we 
also follow the evolution given by the numerical finite-differences 
solution of Eq. 1 140 I I (thick line) and the analytical solution (line 
with empty squares) Eq. I I41I I whose parametres were chosen in a 
such way as to share the same magnetic flux. We used the parame- 



tres /3 = 1.1, t^j-ag/t 



drag I '^weak 



100. The labels in the figure represent 



the instants: (1) t = 0, (2) t = Wt^ea^, (3) t = 
used the same background parametres noc/non = 0.04 



^(drag) _ ^ 
ambip 



We 



process acts in a tendency to spread out the magnetic flux. 
Again, as expected, the total magnetic flux is conserved. 

• At the nuU points, the magnetic field vanishes contin- 
uously, but with a high (possibly inflnite) derivative, as ex- 
pected from the example of the solution given by Eq. (|41|l . 

• Contrary to the weak coupling limit, in Fig. [O] it is 
observed that the magnetic flux is not preserved in each 
one of the regions separated by the magnetic null points. 
Therefore, there is a transfer of magnetic flux through these 
null points, leading to reconnection of magnetic field lines 
even in the absence of Ohmic diffusion, as will be discussed 
further below. 



4-2.3 Fixed singularity 

In order to understand the formation and behaviour of the 
singularities at the null points, we first consider a magnetic 
field with a finite derivative at a null point x = xo, and 
follow the evolution of its derivative p{t) using the Ansatz 

B{x,t)^p{t){x-xo). (44) 

The transport velocity, 

Vc = —^{B^)^-2p\x~xo), (45) 

vanishes at the null point as long as p is finite, so there is no 
reconnection at the null point in this regime. From Eq. (|40p . 
the differential equation for p{t) is dp/dt — 4p'^, with the 
solutionp(t) = p(0)/a /1 — 8p^(0)t, which diverges in a finite 
time too = 1/(80^(0)). [Brandenburg fc Zweibef(|l994l ) found 
that this divergence (which they identified numerically) can 



Figure 8. Evolution according to Eq. 1 140 II of an initial field 
B{x,Q) = exp(-400(x - 0.2)2) ^ 0.7 exp(-100(x - 0.8)^). The 
label numbers show the time progression for different instants, 
(1) i = 0, (2) t = 0.3, (3) t = 3.0 
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Figure 9. Evolution according to Eq. 1 140 I I of an initial field 
B(x, 0) = exp(-400(x - 0.2)2) _ o.7cxp(-100(x - 0.8)2).The la- 
bel numbers show the time progression for different instants, (1) 
t = 0, (2) t = 0.5, (3) t = 45.0 



lead to a stationary solution for the magnetic field. We re- 
derive it by noting that, to have dB/dt = in Eq. (|40p . 
the term in parenthesis on the right-hand side, which cor- 
responds to (minus) the "flux of flux" (amount of magnetic 
flux crossing any point x per unit time), must be uniform in 
space. 



B 



dBf_ 

dx 



'F = constant, 



leading to 



Bix) 



-F{x - xo) 



-,1/3 



(46) 



(47) 



The infinite derivative at x = a;o allows F to remain finite at 
this point, so there is magnetic fiux crossing the singularity 
and causing reconnection, without having included Ohmic 
diffusion in the model. Of course, this simple solution is not 
compatible with our boundary conditions, which were con- 
structed so as to enforce F = at the boundaries. Therefore, 
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for our numerical solutions, the reconnection must produce 
a decrease in the absolute value of the magnetic flux in each 
of the regions separated by the singularity, and therefore in 
the amplitude of the magnetic field. 

In order to explore this behaviour, we will study the 
evolution of an initial magnetic field profile of the form 



B{x, 0) = cos{nx). 



(48) 



This magnetic field profile has a null at a; = 1/2 and is 
an odd function respect to this point, a property that is 
preserved by the evolution according to Eq. (|40|) . Thus, in 
particular, the position of the null point will remain fixed. 
Based on this fact, we try a solution of Eq. (1401) by separation 
of variables in the form 



with f(t) and g(x) satisfying 



and 



d_ 

dx 



idg_ 
dx 



C 



(49) 



(50) 



(51) 



where C is a separation constant. The magnetic flux between 
the left boundary of the system and the null point at x = 1/2 
is 

nl/2 

(52) 



$+(t) = / B{x,t)dx. 
Jo 



By symmetry, the flux in the other half of the interval will 
be 



$"(*)=/ B{x,t)dx = -'i>^{t). 

J 1/2 



If we normalize, 



1/2 



g{x)dx — 1, 



(53) 



(54) 



we wiU have f{t) — <I>^(t). 

In order to solve Eq. (|5ip subject to the boundary condi- 
tions of Eq. (|8}, we use the new variable u{x) — g'^{x), which 
must satisfy the conditions {du/dx)a:=o — {du/dx)x=i = 
and it(l/2) = 0. Thus, Eq. (|5ip is rewritten as 

cPm _ 3 1/3 

dx^ - 2^" ■ 



(55) 



Combining Eq. (|55p with Eq. (|54[) . we obtain the auxiliary 
condition {du/dx)^^i^2 = — 3C/2, which indicates that the 
magnetic field close to the null point has the shape given 
by Eq. (|47p . The numerical solution of Eq. H55p must satisfy 
conditions at different points, which can be intricate from 
the numerical point of view. Thus, we rescale of our variables 
as -u = Uv and x = Xy, where U and X are constants to be 
determined. This allows us to write Eq. (f55]l as 



dy'2 



„l/3 



(56) 



which requires C = (2/3)X-2f/^/^ The boundary condi- 
tions can be set as v{y — 0) — 1, {dv/dy)y^o — 0. From 
the numerical solution of Eq. (|56p we obtain v{yo) = 
at yo = 1.46. Since the null point of u is at xq — 1/2, 
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Figure 10. Evolution of an initial magnetic field B{x, 0) = 
cos(7ra;). The label numbers show the time progression for dif- 
ferent instants, (1) t = 0, (2) t = 0.5, (3) t = 5, (4) t = 50. 
The dashed lines represent the evolution as given by Eq. I|40|l 
while the full line is the asymptotic solution given by Eq. I I49I I. 
It follows that the solution of Eq. 1140 |l converges after an initial 
transient to Eq. 1149 l l. 3'+(t) is the positive magnetic on the left 
of the singularity x = 0.5 [see Eq. | [57t ] and <I>+(0) = I/tt. 



we get X = xo/yo = 0.34. Using u — Uv and the condi- 
tions on the first derivatives, [du/ dx)^^i/2 ~ — 3C/2 and 
{dv/dy)y^yg = po = —1.20 (obtained numerically), we get 
U = ~l/{Xpof = 14.0. Finally, comparing Eqs. (|55p and 
(f56ll we obtain C = (32/3) (yo/p^) = 33.7, which determines 
the full evolution of the magnetic flux in each half of the 
interval, which is given according to Eq. (|50p as 



$+(0) 



1 + 2C ['i>+iO)ft 



(57) 



The flux calculated from the numerical solution of Eq. 1)40^ 
agrees with Eq. (|57p with a percentage error of order 0.4%. 
Eq. (|57l) implies that positive magnetic fiux coming from 
the left of the singularity annhilates with the negative flux 
coming from the right. This is reconnection in the absence 
of Ohmic resistivity. 

In Fig. [To] we observe that the magnetic field as given by 
the numerical solution of Eq. (|40|l converges after a initial 
transient to the field of Eq. (|49p (obtained from the solution 
of Eq. ([55ll. 

This magnetic field profile has a strong divergence of 
the velocity field at the null point, and we are assuming 
that the excess charged particles instantaneously recombine 
when they reach this point, which is not realistic even in the 



"strong coupling limit", in which t^^ak «C tdr 



In other 



words, even if the beta decays can be considered as instanta- 
neous everywhere else in the system, this approximation will 
break down close enough to the singularity. The behaviour 
close to the singularity, including both a fast, but finite re- 
combination rate and a s mall, but finite O hmic diffusivity, 
was discussed in detail bv lHeitsch fc Zweibel (|2003a ibi). 



5 CONCLUSIONS 

We have studied the asymptotic magnetic field evolution 
promoted by ambipolar diffusion in a one-dimensional geom- 
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etry for two opposite, limiting regimes. In the weak coupling 
limit, in which neutral and charged particles drift easily with 
respect to each other, the bottleneck for the evolution is the 
conversion of one species into another, which is required in 
order to eliminate the charged-particle pressure gradients 
caused by the magnetic field, which impede its evolution. 
In the strong coupling limit, conversions are easy, but the 
inter-particle collisions are the corresponding bottleneck. In 
molecula r clouds, the second regime appears to be generally 
relevant (|Shull 19831 : [Brandenburg fc Zweibelll 19941) . Neutron 
stars in their hot, early phase will also be in the strong- 
coupling regime, and evolve to the weak coupling regime as 
they cool. 

In the weak coupling limit, the magnetic field evolution 
is described by a non-linear, partial integro-differential equa- 
tion, while in the strong coupling limit this evolution is de- 
scribed by a non-linear diffusion equation. We made numer- 
ical simulations of the evolution of different initial magnetic 
field profiles in each of these limits and found agreement 
between our numerical results and some analytic solutions 
that can be found for these differential equations. 

From our results, we infer that, in both limits, the am- 
bipolar diffusion process operates in a tendency to spread 
out the magnetic flux, but contrary to the normal Ohmic dif- 
fusion this process asymptotically produces singular points 
with sharp magnetic field gradients. These sharp gradients 
develop around those points where the magnetic field is null, 
and separate regions of magnetic fields with opposite signs. 
We observe some generic properties of this process, as fol- 
lows: 

In the weak coupling limit, the resulting discontinuities 
can be modelled as step solutions [Eqs. (|29|) . (|24p ]. The 
asymptotic magnetic field is spatially uniform in each of 
the regions separated by these singularities, and its abso- 
lute value (and thus the magnetic pressure) is the same in 
each region. Ambipolar diffusion by itself does not cause 
magnetic fiux transfer (and thus reconnection) across the 
singularities, but the associated current sheets will easily be 
dissipated by Ohmic diffusion, so reconnection will occur in 
a realistic system. In the strong coupling limit, at the sin- 
gular points the magnetic field vanishes continuously but 
with infinite spatial derivative. Ambipolar diffusion acts in 
a tendency to spread out the magnetic flux, but, contrary 
to the weak coupling limit, the magnetic flux is not pre- 
served in each one of the regions separated by the magnetic 
null points. Therefore, there is a transfer of magnetic flux 
through these null points, i. e. reconnection without Ohmic 
resistivity (but see lHeitsch fc Zweibe il l2003al lbh. 

The main limitation in applying the present formalism 
to realistic systems (either neutron stars or molecular cloud 
cores) is the very restrictive, one-dimensional geometry. An 
extension to more realistic geometries (i. e., axial symmetry) 
will be attempted in further work. 
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